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Abstract 

We have mapped the molecular-atomic transition in liquid hydrogen using first principles molec- 
ular dynamics. We predict that a molecular phase with short-range orientational order exists at 
pressures above 100 GPa. The presence of this ordering and the structure emerging near the dis- 
sociation transition provide an explanation for the sharpness of the molecular-atomic crossover 
and the concurrent pressure drop at high pressures. Our findings have non-trivial implications 
for simulations of hydrogen; previous equation of state data for the molecular liquid may require 
revision. Arguments for the possibility of a 1'^* order liquid-liquid transition are discussed. 

PACS numbers: 62.50.-p,61.20. Ja,64.70.dj,71.22.+i 
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Constructing an accurate picture of high pressure hydrogen continues to present a chal- 
lenge for both experiment and theory. In recent years, many studies have focused on the 



behavior of hydrogen at elevated temperatures. Topics of interest inc 




flflfl 



ude the turnover of 



a semiconductor 



the melting line [1|, y, [a , the molecular-atomic transition 
to metallic transition |9|, [lO| , as well as a host of temperature-induced dissociation and ion- 
ization phenomena. Despite considerable attention, many important questions about the 
chemical and physical properties of hydrogen remain unanswered. A complete description 
of the molecular-atomic transition has not been established. The relation between metalliza- 
tion and dissociation is still not fully understood. Effects due to impurities (e.g. helium) [3] 
and the accuracy of equation of state (EOS) data are particularly important for determining 



the structure of gas giants such as Jupiter 



ill. 



This Letter is focused on a topic which has hitherto received little attention - the structure 
of the dense hydrogen liquid above the melting curve and near the dissociation transition. 
Based on first principles molecular dynamics (FPMD) simulations, we predict a liquid phase 
that exhibits complex short-range orientational order. This finding is significant, as the 
properties of this phase provide a physical explanation for some of the phenomena mentioned 
above. 



We have performed FPMD simulations using DFT-GGA |16l.ll7l|. mapping the molecular- 
atomic transition over a large pressure (P) and temperature (T) range (Fig. [1]). Simulations 
were carried out using the Born Oppenheimer implementation in the CPMD code [isl]. In 
this method, the electronic density is optimized to its ground state at each ionic step. We 
used a local TrouUier- Mart ins pseudopotential with a 100 Ry planewave cutoff and P-point 
sampling of the Brillouin zone. Molecular dynamics (MD) simulations were carried out in the 
NVT ensemble, where a Nose-Hoover thermostat was used to control the ionic temperature. 
The supercell in our simulations consisted of 256 atoms. We performed checks on larger cells 
of 512, 768, and 1024 atoms; discussion of these results follow. An MD time-step of 16 a.u. 
(0.387 fs) was used for all simulations, with checks performed using a time-step of 8 a.u. 
Observables such as the pair correlation function, stress tensor, self-diffusion coefficient, and 
spacial distribution functions (discussed below), were all found to be well converged. In all 
cases, simulations were allowed to run for at least 2 ps and some as long as 15 ps. For the 
densities considered here, the system rapidly reaches equilibrium and these simulation times 
are sufficient. Calculations of the polarizability tensor were produced by extracting a large 
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FIG. 1: Phase diagram of dense hydrogen. Triangles are upper and lower boundaries of the 
computed molecular-atomic transition. The T range over which the transition takes place decreases 
with pressure, consistent with previous observations [^,0, ^, 3]. Experimental 0, Q] and 
theoretical data (extrapolation for P > 200 GPa) of the melting line are shown along with 
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the semiconducting-metallic transition 
estimated to be 30% pji]. QMC simulations 
and 364 K. 



. . Uncertainties in the T of this measurement were 
I5I confirm the stability of the liquid near 320 GPa 



number of independent configurations from equilibrated molecular dynamics trajectories and 
using linear response to calculate the polarizability of each. 

To map the liquid phase diagram (Fig. [1]), we calculate the survival probability, n(r), of 
paired atoms for different P, T conditions (over 150 samples, each corresponding to a fully 
equilibrated molecular dynamics trajectory). A pair is defined to be two H atoms which 
are mutual nearest neighbors. Given a pair at time t = 0, we calculate the probability for 
the same pair to exist at t = r (see Supplimentary Fig. 1). For a purely molecular system, 
n(r) = 100%, while for a purely atomic one II{t) — > 0% (see Supplimentary Fig. 2). We 
found that our results are qualitatively insensitive to the parameter r (we propose a standard 
of 10 H2 oscillations). This criterion relates directly to the rate of dissociation in the liquid, 
and naturally accounts for bond breaking and reforming characteristic of high pressure 
liquids. The molecular-atomic transition based on II{t) = 50% is shown as a function of 
P in Fig. [Hand as a function of the density parameter (defined by V/N = |7r(rsao)^5 
where V is the volume, N the number of electrons, and is the Bohr radius) in Fig. [2](b). 
Linear extrapolation of the transition line suggests that ~ 1.25 is the maximum density 
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where molecules would be stable in a T = K liquid (the existence of which has been 
previously proposed |l, [lol, Q). The inclusion of quantum zero-point effects is expected 
to shift the onset of dissociation to lower densities. On the other hand, the stability of 
hydrogen pairs in a solid can be enhanced due the presence of intermediate-range order and 
Friedel oscillations 21|. 



To relate the structural changes to the dielectric properties of the hquid, we have com- 
puted the static dielectric constant (average of the diagonal elements of the polarizability 
tensor). As shown in Fig. [21 it increases by several orders of magnitude when the system 
is brought through the dissociation transition. As the probability of molecular survival 
decreases, the liquid becomes increasingly polarizable. This change is consistent with a 
transition from an electronically insulating to a conducting state - a metal or a thermally 
activated semiconductor. 
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FIG. 2: (a) Static dielectric constant of the liquid calculated along compression isotherms. Error 
bars indicate standard deviation of results from different atomic configurations. The pressure scale 
corresponds to the 1500 K isotherm. Significant changes in the dielectric properties of the liquid 
begin at the density corresponding to the onset of dissociation (b). Definition of symbols is the 
same used in Fig. [H 



The appearance of a maximum in the melting line of hydrogen near 82 GPa has been 
explained by a softening of intermolecular interactions in the liquid similar to what 
was observed in the orientationally-ordered phase of the solid (phase III) j22|. This has 
prompted us to investigate the structure of the molecular liquid. We have discovered a new 
region in the molecular fluid in the vicinity of the melt line turnover that exhibits short 
range orientational ordering. It is revealed in the spatial distribution function, SDF{r,6). 
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Given an H2 molecule, SDF{r, 9) is defined as the probability for finding a hydrogen atom 
at a distance r from the molecular center of mass and at an angle 9 with the molecular axis, 
normalized by the average H number density. 

In Fig. [31 we show a graphic of the SDF for two densities along the same 1000 K isotherm; 
for clarity only the contributions from nearest-neighbors are shown. In the low density case 
(vertical plane), a molecule's nearest neighbor has an equal probability of being found at 
any 9. At high pressure (horizontal plane), there is an orientational correlation between 
neighboring molecules. Particles are less likely to be found at the molecular poles, and more 
likely to be found near the equators. This can be further seen in the inset of Fig. [3l where we 



have integrated over the radial degree of freedom, JSDF{r,9) dr = SDF{9). To describe 
the evolution of the structural transition under compression, we define an order parameter 
as a = 'SDF{9 = 90) fSDF {9 = 0). A plot of a as a function of r, along the 1000 K 
isotherm (Fig. H]) indicates a rapid increase in the vicinity of the previously reported l|, [2I, 
maximum in the melting curve. It is this change in the structure of the liquid, driven by a 
change in the intermolecular interactions, that is responsible for the turnover in the melting 
line. 
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FIG. 3: Structural order parameter, a, and spatial distributions functions, SDF{r, 9) and SDF{9) 
(see text for definitions), along the 1000 K isotherm in molecular H. The density where the melting 
curve has a maximum is indicated on the plot; it is in the region where a begins to increase rapidly. 
The graphic shows nearest neighbor contributions to SDF{r, 9) at low (r^ = 2.40, P = 6 GPa) 



and high (r^ = 1.45, P = 170 GPa) density. The SDF{0) further illustrate changes in the angular 
distribution with density. 



A detailed quantitative description of the emerging orientational order will be presented 
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in a follow up paper. Such analysis is expected to have relevance for identifying the finite 
(and possible zero) T phases of the solid; similarities with phases II and III of H are plausible. 
Here we mention only that in addition to the described dependence on 9, there is a very 
strong tendency for neighboring molecules to lie in plane. For 9 ~ 90°, the angle between 
molecular axes tends to be about 15° (for different 9 it depends slightly on whether the 
molecules are tilted towards or away from each other). 

One of the most prominent features of dense H is the sharp transition from molecular to 
atomic liquid. At sufficiently high density, the onset of dissociation leads to an EOS with a 
negative slope, i.e. dP / dT\v=const < 0. There has been considerable debate as to whether 
the abrupt change in P over a narrow T range is due to a discontinuous (P* order) transition 
in the liquid and how it is affected by the presence of impurities (e.g. He). In what follows, 
we show that the anomalous features in the EOS, its density and impurity dependence, and 
its sensitivity to computational parameters can be understood on the basis of the changes 
in the liquid structure across the dissociation transition. 

In order to describe the liquid structure in the presence of both atoms and molecules, 
we decompose SDF{r, 9) into atomic and molecular components. In both cases, SDF{r, 9) 
is calculated with respect to reference molecules center of mass. For the atomic SDF{r,9), 
statistics for r and 9 are collected for single atoms only, while for the molecular one, the atoms 
at r and 9 are paired. In Fig. H] we show such a decomposition at conditions corresponding 
to a large degree of local angular structure {a ~ 10). At these P and T, the liquid has 
just begun to dissociate (H(r) ^ 99%). The molecular SDF is similar to that of a purely 
molecular liquid, with a peak near the reference molecule's equator. However, the atomic 
SDF shows that unpaired atoms tend to shift towards the molecular poles. The atomic 
^DF{9) is peaked at ^ ^ 50° (and 130°) versus 90° for the molecular JDF{9) [Fig. H^b)]. 
Similar analysis at low density shows no such angular dependence. 

At low P, where there is no short-range orientational order, the H2 molecules are freely 
rotating and their packing is similar to that of spheres. Upon dissociation, the problem 
becomes that of spheres of different sizes, however, there is no apparent optimization of 
packing. On the other hand, because of covalent bonding, the effective volume of H2 is less 
than that of two isolated hydrogen atoms. For these reasons, the dissociation transition at 
low density is characterized by dP/dT\v=const > 0. 

At high P, the molecules are no longer freely rotating and the packing problem becomes 
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FIG. 4: Molecular and atomic spatial distribution functions, SDF{r,6) (a) and SDF{6) (b); see 
text for definitions. For clarity, only contributions from nearest neighbor atoms are included. 
Atoms forming molecules are more likely to reside near the reference molecule's equator. The 
distribution of unpaired atoms is more even and indicates that they are most likely to be found at 
e » 50° and 9 « 130°. 

that of prolate spheroids (see graphic in Fig. |4]). As seen Fig. HI the molecular orientational 
order creates voids between them, which can be filled in by single atoms after dissociation. 
There will be a drop P upon dissociation when the optimization in packing is sufficient to 
compensate for the higher effective volume of single atoms compared to that of molecules. 



T 



lis picture is consistent with the fact that the sharpness in the EOS increases with density 



12l | - our order parameter follows the same trend. Furthermore, it has been noted 
that diluting the system with neutral impurities such as helium 7| softens the transition. 
We suggest that performing SDF analysis on high density mixtures of molecular hydrogen 
and helium (or other noble gases) will reveal similar angular local order, with the dopant 
atoms residing within the voids of the liquid. Finally, we note that in a one-dimensional 
hydrogen liquid, where angular order is not possible, dP/dT\v=const must always be positive. 

Packing considerations also provide an explanation for the large sensitivity of the tran- 
sition to the number of atoms, A^, in the simulation supercell. We find that the degree of 
dissociation and related properties depend on in a non-intuitive way (similar observations 
have been made by Morales jisl). At high densities, II{t) oscillates as A^ is varied from 256, 
512, 768, 1024. Surprisingly, simulations conducted with 256 atoms are in better agreement 
with those done with A^ = 1024 than with A^ = 512 (n^=i28 = 0.3338, n^=256 = 0.98 42, 
n;v=5i2 = 0.8328, n7v=768 = 0.9660, and n7v=io24 = 0.9656 for r, = 1.45, T = lOOOif). This 
effect cannot be removed with a finer fc-point sampling of the Brillouin zone. With A^ = 128 
the error in P is as large as 8% at rs=1.40, T=750 K. 
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The origin of this initially puzzling behavior is that at a given density, defines the 
physical dimensions of the simulation box. Depending on how these dimensions relate to 
the correlation length of the liquid, different sized cells will have a slight bias towards one 
phase or another. Ideally one should use a value large enough that this effect is no longer an 
issue. Given the current cost of using large cells, we report values corresponding to = 256, 
as they are in relatively good agreement with N = 1024 (the oscillation in n(r) is quite 
damped in this larger cell). Simulations performed with = 256 tend to favor the molecular 
phase relative to that of = 1024; this should partially correct for the underestimate of 
the dissociation barrier due to approximations of the DFT exchange-correlation potential. 

The existence of short range angular structure in the liquid, specifically the different 
arrangement of molecules and atoms, suggests that something similar might also occur in 
the solid phase. In future searches for solid state structures it would be worthwhile to include 
variations of a two-component solid, comprised of molecules and atoms. Indeed, work by 



Pickard and Needs [23=] indicates that mixed layered structures, comprised of molecules and 
atomic graphene-like sheets, remain energetically competitive over a wide range of pressures. 

We now turn our discussion to the possiblity of a discontinous {i.e. 1^* order) transition 
in the liquid. Based on our sampling of the phase diagram, the molecular-atomic transition 
appears to be continuous up to 200 GPa. We do not find a fiat region in the P{V) EOS - a 
signature of a order transition. Our data imply one of three possibilities: the transition 
is not 1** order, it takes place outside of the P, T space we considered, or it occurs over 
a range of densities too small to be resolved by our sampling. We can use our results, 
however, to establish a connection between the microscopic and macroscopic properties 
characterizing the transition and provide insight for the physical mechanism that could lead 
to it being 1** order. Let 6V be the reduction in volume (due to packing) that can be 
achieved by dissociating a single molecule, while keeping P and T constant (the resulting 
state need not be in equilibrium). If a 1'^* phase transition exists at these P and T, we 
must have —PAV = AU — TAS. Here AV^, AU and AS are the volume, energy and 
entropy differences across the transition. If it results in the dissociation of molecules, 
then we define Ed = AU/rid, which has the physical meaning of dissociation energy. A 1** 
order transition will take place if UdSV = —AV, which means 6V ~ E^/P (assuming that 
TAS/P can be neglected, especially at high P and low T). We note that E^ decreases under 
compression due to screening, while SV increases due to more pronounced orientational 
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order. Thus, convergence of these two terms, SV and Ed/P, is increasingly hkely at high 
P and low T. Conclusive determination of the existence of a critical point will require very 
fine sampling of the phase diagram below 1000 K and above 200 GPa. 

We have mapped the molecular-atomic transition of liquid hydrogen over a large pressure 
range. Our transition line correlates well with changes that can be observed in the electronic 
properties of the system such as the static dielectric constant. We predict that the liquid 
demonstrates significant short range orientational ordering. Its development coincides with 
the turnover in the melt line. The existence of this structural order is responsible for the 
large dissociation-induced P drop in the EOS. Furthermore, it provides an explanation for 
the significant finite size effects that are present in this system. A 1*** order transition in the 
liquid, if it does exist, will likely occur at T < 1000 K and P > 200 GPa. 
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